EHT Imaging Tutorial 


Andrew Chael & Katie Bouman 


Imaging with the Event Horizon Telescope 


Measurements 





Infinite Number 
of Possibilities 





Imaging Methods 


SQUEEZE, BSMEM, Closure-only, Sparse Imaging, CHIRP, Closure-only 





eht-imaging Python Library 
github.com/achael/eht-imaging 
* Generate Data 


* Plot Data/Results 


* Closure-only Imaging 
* Time-Variable Imaging 


ә Scattering Mitigation 


Step 1: Exploring the VLBI 
Imaging Website 


VLBI Imaging Website  vlbiimaging.csail.mit.edu 


2” VLBI Reconstruction Dataset 


A Dataset Designed to Train and Test Very Long Baseline Interferometry Image Reconstruction Algorithms 





Welcome to the VLBI Reconstruction Dataset! 


The goal of this website is to provide a testbed for developing new VLBI reconstruction algorithms. By supplying a large set of easy to understand training and testing data, we 
hope to make the problem more accessible to those less familiar with the VLBI field. Specifically, this website contains a: 


e Large set of synthetic training data for many different VLBI arrays and targets 
e Set of real data measurements provided in the same standard format 


e Standardized data set for testing VLBI Image Reconstruction Algorithms 

e Online quantitative evaluation of algorithm performance on simulated testing data 

e Qualitative comparison of algorithm performance on the reconstruction of real data 

e Online form to easily simulate realistic data using your own image and telescope parameters 


VLBI Imaging Website  vlbiimaging.csail.mit.edu 


Standardized dataset of real & synthetic data 


Over 5000 synthetic measurement sets: 
14 Array Configurations, 96 Source Images, 4 Noise Levels 


VLBI Imaging Website  vlbiimaging.csail.mit.edu 


Automatic Quantitative and Qualitative Evaluation 


VLBI Imaging Website  vlbiimaging.csail.mit.edu 


Online form to easily simulate realistic data using 
user-specified parameters 


Step 1: Generating Data on the 
VLBI Imaging Website 


Selecting/Uploading an Image 


Step 1: Select Image of the Emission 


Sorry about this! 
We will fix the 


Ninconsistency 
Select or upload an image that you would like to observe and specify the total flux density of the emission. ve ry SOON 


Total Flux Density (Janskys) 


CLICK TO 
UPLOAD YOUR 
DWN GRAYSCALE 
PNG/JPEG IMAGE 
under 512 x 512 pixels 


> 
Add Image To Dataset sgr " model 
When using these images cite 


(Broderick et al., July 2011) 





Spin: 096 
Inclination: 89 ° 


Ҹә 4» 


Rotation (Degrees): 








Natural User Uploaded 


Celestial 


Choose Another Choose Another Choose Another 
Image Image Image 


Selecting Target Location and Field of View 


Step 2: Select Direction and FOV 


Identify the direction to the target source. Right ascension should be in the form HH:MM:SS.SS for hours, minutes, and seconds and declination should be in the form 
DD:MM:SS.SS for degrees, arcminutes, and arcseconds. Field of view is specified in arcseconds. Warning: You must choose coordinates such that your region will be observable 
from your observatory site (the first telescope you specify below) at the start time that you specify, otherwise the resulting output will be incorrect. 


Field Of View Center: Right Ascension (HH:MM:SS.ss)| — 174540041 | Declination(DD:MM:SS.SS) © 290028118 
Field Of View Size: Right Ascension (areseconds) ^ 000016. | Declination (arcseconds) — 000016 d 








Selecting Telescopes 


Step 3: Specify Telescope Array 


Add the telescope locations and intrinsic parameters that you would like to use to simulate data 


Initilization: Select a pre-loaded telescope 


Name: Unique name for each telescope station (up to 12 characters) 


East Longitude/Latitude: East longitude and latitude of the array center. For locations less than 180 degrees west of Greenwich a Minus sign should precede the longitude entry. 


X/Y/Z Position: Absolute X, Y, Z coordinates of each station (in meters) relative to the center of the Earth 
Lower/Upper Elevation: Lower and upper elevation limits of the of the antenna in degrees 
SEFD: System equivalent flux denisty of the antenna 


Diameter: Antenna diameter in meters 


Initilization Name East Longitude Latitude X-Position 

ALMA + ALMA -67:45:11.4 -23:01:09.4 2225037.1851 
SMT SMT -109:52:19 32:42:06 -1828796.2 

m LMT v LMT -97:18:53 18:59:06 -768713.9637 
SMA v SMA -155:28:40.7 19:49:27.4 -5464523.4 
PV E PV -3:23:33.8 37:03:58.2 5088967.9 
PDB $ PDB 05:54:28.5 44:38:02.0 4523998.4 
SPT v SPT -000:00:00.0 -90:00:00 0 


ADD TELESCOPE DELETE SELECTED 


Y-Position 


-5441199.162 


-5054406.8 


-5988541.7982 


-2493147.08 


-301681.6 


468045.24 


Z-Position 


-2479303.4629 


3427865.2 


2063275.9472 


2150611.75 


3825015.8 


4460309.76 


-6359587.3 


Data Collection Settings 


Step 4: Specify Date and Time Data is Collected 


Specify the time of when you would like measurments to be taken, and the time interval between measurements. 


Start Time: Specify the time of your first observation in Universal Time (UT). The required format is "YYYY:ddd:hh:mm:ss" where YYYY is the year, ddd is the day number (e.g., December 31 is 
day 365); hh is the UT hour, mm is the UT minute, and ss is the UT second. 

Scan Duration: The length of a continuous scan in seconds 

Interval Length: The time in seconds between successive scans 


Number of Samples: The number of successive scans of this type 


Start Time (UT) Scan Duration (seconds) Interval Length (seconds) Number of Samples 


ADD DATA COLLECTION DELETE SELECTED 


Data Collection Settings 


Step 5: Specify Collection Parameters 


Specify the center frequency and width of the observing channel in MHz. 


Center Frequency (MHZ! 297297 | Bandwidth (MHz) дәв 


Specify your integration time in seconds (sometimes referred to as "dump time" or "record length"). This is not the total duration of your observation, but rather the sampling and 
recording interval of the data. 


Integration Time (seconds) ор | 


What Kinds of Noise Can We Add? 


Thermal Noise "Atmospheric Phase Error Systematic Gain Error 


Selecting Types of Noise Added 


Let's JUST add Thermal Noise 


Step 6: Add Noise and Generate Data 
Simulate Without ANY Noise 
Уј) Simulate Without Atmospheric Phase Errors 


«| Simulate Without 5 % Gain Error 


And now we generally generate our data.... 


But if so many people submit at the same 
time we will probably bog down the 
machine.... 


So for now, please download pre-computed 
data and later you can generate it yourself 


vibiimaging.csail.mit.edu/myDataResults 6312 


Click Here to Download Data 


10 1e10 


UV Coverage 


GOG 








xx SMT-ALMA 
xx LMT-ALMA 
xx LMT-SMT 
xx SMA-ALMA 
xx SMA-SMT 
xx SMA-LMT 


0.5 


xxx PV-ALMA 
xx PV-LMT 
xx PDB-ALMA 
xx PDB-PV 


0.0 


x 
xXx SPT-ALMA 
xx SPT-SMT 
xx SPT-LMT 
xx SPT-SMA 
xx SPT.PV 
xx SPT-PDB 





V (Wavelengths) 


-0.5 





U (Wavelengths) ue 


Click Here to View the Telescope and Target Source Parameters 





Visibility Amplitudes Visibility Phase Plotted on Real Imaginary Plane 
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What does the downloaded zip file provide? 


uu Information to Reproduce Data 
url.txt 
v EE data 


“| sgraimage.mat 


sgraimage.oifits | 
sgraimage.txt ata in a number of formats 
saraimage.uvfits 
UVdata.txt 
Y | | plots 
— sgraimage amp MAPS.png 
— Sgralmage, phase. MAPS.png Plots to help you understand the data 
| sgraimage visi MAPS.png 
|| uvcov.png 
“| sgraimage.fits 
B sgraimage.png 





iginal Images in FITS and PNG 


What are these data formats? 


What are the data formats and how do I use them? 


We use OIFITS, MAT, and ASCII data formats. We describe OIFITS and MAT below in detail. 


OIFITS 
The primary data format that we have chosen to use is OIFITS. OIFITS is a standard for exchanging data for Optical (Visible/IR) Interferometry, and is based on the FITS 


Standard. Since mm/sub-mm VLBI shares a lot of similarities to optical interferometry, this format is better suited for mm/sub-mm measurements than UVFITS. More information 
on the OIFITS format can be found here. We list the variables described in this paper in the tables below. 


We provide a number of tools that may be useful in reading and writting in the OIFITS format: 


e Paul Boley has written a OIFITS Python module that you can download here 
e Python code that can be used to write an OIFITS file from an output MAPS text file can be downloaded here 
e Python code by Andrew Chael to extract information from OIFITS and write it to a text file can be downloaded here 


OI T3 Variables Description Units 
T3AMP Triple-product/Bispectrum amplitude Jansky? 
T3PHI Triple-product phase Degrees 
T3AMPERR Standard deviation of error in triple product amplitude Jansky? 
T3PHIERR Standard deviation of error in phase Degrees 
U1COORD u coordinate of baseline AB of the triangle meters 
U2COORD u coordinate of baseline BC of the triangle meters 
V1COORD v coordinate of baseline AB of the triangle meters 
V2COORD v coordinate of baseline BC of the triangle meters 
STA INDEX Station numbers contributing to the data 

INT TIME Integration time seconds 


MJD Modified Julian Date 


Step 3: Loading and Inspecting 
Data 


In an ipython window: 


import numpy as np 


import ehtim as eh 


Load the observation file we generated 


obs = eh.obsdata.load uvfits('./data/sgraimage.uvfits') 


Can also load custom text format, 
oifits, and MAPS output 


Look at plots! UV coverage 


obs.plotall('u','v', conj=True) 


Shows both u,v and -u,-v 





Look at plots: Visibility amplitudes 
Other possible fields include 
obs.plotall("uvdıst”, ‘amp’ je Snr”, "sigma", “qamp”, 


// uamp // : // vamp // : Ы diki 





0.6 
uvdist lelü 


Look at plots: Baseline phase over time 


obs.plot biI('SMA','ALMA', "phase") 


SMA - ALMA 





10 
GMST (h) 


Look at plots: Closure phase over time 


obs.plot cphase('LMT', 'SPT', 'ALMA') 





GMT (h) 


Take a look at the dirty beam and clean beam 


o ARGET NAME DB MJD О 227.30 GHz 


Image Parameters = 


FT of the sparse u-v coverage = 


npix = 128 
fov = 200*vb.RADPERUAS 
Dirty Beam 


— 





dbeam = obs.dirtybeam(npix, fov) = 
db eam. cl 1 S p | a V ( ) "BU uM iu SEN Mn EEUU “əl -100.0 


Gaussian fit to the beam component 


00 og ARGET NAME MJD О 227.30 GHz 

0.0054 
75.0 

0 





Clean Beam 
cbeam = obs.cleanbeam(npıx,fov) 
cbeam.display() 


0.0048 





Take a look at the dirty image 


dim = obs.dirtyimage(npix, fov) 
dim.display() 
TARGET NAME MJDO 227.30 GHz 





0.00075 


0.00060 


Sky Image convolved with Dirty Beam 222 


10.00045 
i." VW 25.0 
ER ä 0.00030 
Q 
- 0.0 
" 0.00015 
2 
© 
@ -25.0 


0.00000 


—0.00015 


—0.00030 





-100,0 —— — "ya 
100.0 75.0 50.0 25.0 0.0 -25.0 -50.0 -75.0 -100.0 


Relative RA (pas) 


What is the array resolution? 


Clean Beam Gaussian parameters 
= (FWHM max, FWHM min, Position angle) 
beamparams = obs.fit beam() 
Save for use in restoring beam convolution 
"Maximum" resolution 
res = obs.res() 





1/longest baseline — can use in circular 
restoring beam 


Step 4: Produce an Image 


Generate a prior image 


Image Parameters 
npix = 128 
fov = 200*vb.RADPERUAS 


Gaussian Prior 


zol = 2.5 ee O The zero baseline flux 
prior fwhm = 100*eh.RADPERUAS «4— FWHM of our circular Gaussian Prior 


gaussparams = (prior fwhm, prior fwhm, 0.0) 
emptyprior = eh.image.make square (obs, npix, fov) 
gaussprlor = empiyprlor.add Gauss (201, gaussparams) 


gaussprior.display () 


Use MEM with complex visibilities 


out = eh.imager əb”. 


Initial Image 


Prior Image Total flux constraint 


нас, gaussprior, 201, 
“TAS % elena dla), 


sl="gs",maxit=100) 


Pa step: 101 x?: 7.129543 


Data Term & 


Weight 


Relative Dec (pas) 


-100.0 


100.0 H of iterations 


75.0 


50.0 
egularizer Type 

other options: 
tv”, “11”... 





100.0 75.0 50.0 25.0 0.0 -25.0 -50.0 -75.0 -100.0 
Relative RA (uas) 


Blur and restart 


M Fractional beam size 


outblur = out.blur gauss(beamparams, 0.2) 

We decreased data weight to prevent ovei 
out = outblur 
out = eh.imager func(obs, out, out,zbl, dl-"vis", alpha dl-10 , 


sl="gs", maxit=150) 


step: 151 x” : 0.985278 


s) 


Relative Dec (ya 
o 
о 





-100,0 
100.0 75.0 50.0 25.0 0.0 -25.0 -50.0 -75.0 -100.0 


Relative RA (yas) 


Final images — save to file 


| “ 


Final “restored” image 


OlEED Lu” = ONE Blur çi ek Bom 0,35) 


Display results 


out.display () 
outblur.display() 


Save to FITS 


imageout.save fits("./sgraim.fits") 
outblur.save fits('./sgraim blur.fits') 


Relative Dec (uas) 


Relative Dec (gas) 


-100.0 
100.0 75.0 50.0 25.0 0.0 -25.0 -50.0 -75.0 -100.0 


Relative RA (yas) 





-100.0 
100.0 75.0 50.0 25.0 0.0 -25.0 -50.0 -75.0 -100.0 


Relative RA (jas) 








0.0032 


0.0028 


0.0024 


0.0020 


/pixel 


0.0016 3 


0.0012 


0.0008 


0.0004 


0.0000 


0.0027 


0.0024 


0.0021 


0.0018 


0.0015 


/pixel 


0.0012 = 


0.0009 


0.0006 


0.0003 


0.0000 


Look at fit to data - Amplitudes 


eh.plotall obs im compare(obs, out, "uvdist", "amp") 


comp plots.py 

has functions to overplot 
data from different 
observations 





1.0 
uvdist lelO 


Look at fit to data - Phases 


eh.plotall obs im compare (obs, out, "uvdist", "phase") 





Step 5: Generating Data with 
Atmospheric Noise 


Selecting/Uploading an Image 


Step 1: Select Image of the Emission 


Select or upload an image that you would like to observe and specify the total flux density of the emission. 


Total Flux Density (Janskys) Rotation (Degrees):| 180) 


CLICK TO 
UPLOAD YOUR 
PNG/JPEG IMAGE 
under 512 x 512 pixels 


> 
Add Image To Dataset sgr " model 
When using these images cite 


(Broderick et al., July 2011) 
Spin: 076 Choose Another Choose Another Choose Another 
Image Image Image 


Inclination: 89 ° 





Natural User Uploaded 


Celestial 


Ҹә 4» 


Selecting Target Location and Field of View 


Step 2: Select Direction and FOV 


Identify the direction to the target source. Right ascension should be in the form HH:MM:SS.SS for hours, minutes, and seconds and declination should be in the form 
DD:MM:SS.SS for degrees, arcminutes, and arcseconds. Field of view is specified in arcseconds. Warning: You must choose coordinates such that your region will be observable 
from your observatory site (the first telescope you specify below) at the start time that you specify, otherwise the resulting output will be incorrect. 


Field Of View Center: Right Ascension (HP:MM:SS.SS) 12:30:49.423382 | Declination(DD:MM:SS.SS) (12232804366 — 
Field Of View Size: Right Ascension (areseconds) ^ 000016. | Declination (arcseconds) — 00016 d 








Selecting Telescopes 


Step 3: Specify Telescope Array 


Add the telescope locations and intrinsic parameters that you would like to use to simulate data 


Initilization: Select a pre-loaded telescope 


Name: Unique name for each telescope station (up to 12 characters) 


East Longitude/Latitude: East longitude and latitude of the array center. For locations less than 180 degrees west of Greenwich a Minus sign should precede the longitude entry. 


X/Y/Z Position: Absolute X, Y, Z coordinates of each station (in meters) relative to the center of the Earth 
Lower/Upper Elevation: Lower and upper elevation limits of the of the antenna in degrees 
SEFD: System equivalent flux denisty of the antenna 


Diameter: Antenna diameter in meters 


Initilization Name East Longitude Latitude X-Position 

ALMA + ALMA -67:45:11.4 -23:01:09.4 2225037.1851 
SMT SMT -109:52:19 32:42:06 -1828796.2 

m LMT v LMT -97:18:53 18:59:06 -768713.9637 
SMA v SMA -155:28:40.7 19:49:27.4 -5464523.4 
PV E PV -3:23:33.8 37:03:58.2 5088967.9 
PDB $ PDB 05:54:28.5 44:38:02.0 4523998.4 
SPT v SPT -000:00:00.0 -90:00:00 0 


ADD TELESCOPE DELETE SELECTED 


Y-Position 


-5441199.162 


-5054406.8 


-5988541.7982 


-2493147.08 


-301681.6 


468045.24 


Z-Position 


-2479303.4629 


3427865.2 


2063275.9472 


2150611.75 


3825015.8 


4460309.76 


-6359587.3 


Data Collection Settings 


Step 4: Specify Date and Time Data is Collected 


Specify the time of when you would like measurments to be taken, and the time interval between measurements. 


Start Time: Specify the time of your first observation in Universal Time (UT). The required format is "YYYY:ddd:hh:mm:ss" where YYYY is the year, ddd is the day number (e.g., December 31 is 
day 365); hh is the UT hour, mm is the UT minute, and ss is the UT second. 

Scan Duration: The length of a continuous scan in seconds 

Interval Length: The time in seconds between successive scans 


Number of Samples: The number of successive scans of this type 


Start Time (UT) Scan Duration (seconds) Interval Length (seconds) Number of Samples 


ADD DATA COLLECTION DELETE SELECTED 


Data Collection Settings 


Step 5: Specify Collection Parameters 


Specify the center frequency and width of the observing channel in MHz. 


Center Frequency (MHZ! 297297 | Bandwidth (MHz) дәв 


Specify your integration time in seconds (sometimes referred to as "dump time" or "record length"). This is not the total duration of your observation, but rather the sampling and 
recording interval of the data. 


Integration Time (seconds) ор | 


Selecting Types of Noise Added 


Let's add Thermal & Atmospheric Noise 
Step 6: Add Noise and Generate Data 
Simulate Without ANY Noise 
Simulate Without Atmospheric Phase Errors 


+) Simulate Without 5 % Gain Error 


Janskys 





Visibility Amplitudes 


2 





V (Wavelengths) 


1e9 


UV Coverage 





-0.5 0.0 0.5 


U (Wavelengths) 


Click Here to View the Telescope and Target Source Parameters 





Degrees 


—100 


-150 





Visibility Phase 


3 4 5 6 
Wavelengths 


7 


i 





9 
1e9 





vlbiimaging.csail.mit.edu/myDataResults 3593 


Click Here to Download Data 


xx SMT-ALMA 
xx LMT-ALMA 
*» LMT-SMT 
*x SMA-ALMA 
xx SMA-SMT 
xx SMA-LMT 
xx PV-ALMA 
xx PV-SMT 
xx PV-LMT 
xx PV-SMA 
xx PDB-ALMA 
xx PDB-SMT 
xx PDB-LMT 


xx PDB-PV 





1.0 
1e10 


Plotted on Real Imaginary Plane 


Imaginary 





Step 7: Image with Closure Phase 


Look at the phase errors 


Load the data 


obs = eh.obsdata.load uvfits('./data/m87image.uvfits 


phase 


Baseline Phases 
obs.plotall('uvdist','phase!') 





uvdist 1e10 


No phase, No dirty image! „TARGET NAME МЈО О 227.30 GHz 
(Can't CLEAN without Self-Cal) 


Dirty Image - 
npix = 128 

fov = 150*vb.RADPERUAS 

dim = obs.dirtyimage(npix, fov) 
dim.displayvy() 


0.000032 


0.000024 


0.000016 


0.000008 


0.000000 


Relative Dec (pas) 
o 
o 


—0.00000 


—0.00001 


—0.00002 





-75.0 
75.0 56.2 37.5 188 0.0 -18.8 -37.5 -56.2 -75.0 


Relative RA (jas) 


Closure Phase is preserved 


obs.plot cphase(“SMA", "SMT", 'ALMA') 


SMA - SMT - ALMA 


Closure Phase (deg) 





Array resolution and prior image 


Array Resolution 


beamparams = obs.fit beam() 
res — obs.res() 


Prior Parameters 
npix = 128 
fov = 150*eh.RADPERUAS 


zol = 1.0 In, The zero baseline flux 
prior fwhm = 100*vb.RADPERUAS === FWHM of our circular Gaussian Prior 
gaussparams = (prior fwhm, prior fwhm, 0.0) 


Create the Gaussian Prior 


emptyprior = eh.image.make square (obs, npix, fov) 
gaussprior = emptyprior.add gauss(zbl, gaussparams) 


Image with amplitude and closure phase 


DUC = çi o Ee (Obs, ne, Gi poe, ZOL CUP ux CL; 
alpha Gi=100, alpha d2-50, S1-"gs". mqaxrit-l00) 


step: 101 x^: 20.293294 


From experience, closure 
phase fits faster so we 
decrease its weight 


Relative Dec (pas) 





795,0 56.2 37.5 18.8 0.0 -18.8 -37.5 -56.2 -/5.0 
Relative RA (jas) 


Blur and re-image 


outblur = out.blur gauss(beamparams, 0.5) 
out-outblur 
out = eh.imager func.maxen amp cphase(obs, out, out, zbl, dl-"amp", 


AZ Che — alpha Ge ŞU, alpha 2-45. Mekik, St) 


step: 151 x”: 1.236238 


Relative Dec (pas) 





"75.0 56.2 37.5 18.8 0.0 -18.8 -37.5 -56.2 -75.0 
Relative RA (yas) 


Final images 


| “ 


Final "restored" image 


outblur = out.blur gauss(beamparams, 0.5) 


Relative Dec (jas) 





Display results 450 562 375 188 00 188 .375 m Sena 


Relative RA (pas) 


out.display() 


TARGET NAME MJD О 227.30 GHz 
outblur.display() 75.0 J 





Save to FITS 
imageout.save fits('./M87im.fits') 
outblur.save fits('./M87im blur.fits') 


0.0018 


0.0015 


0.0012 


Relative Dec (pas) 


0.0009 


0.0006 


0.0003 





-75.0 
75.0 56.2 375 188 0.0 -188 -37.5 -56.2 -75:0 


Relative RA (yas) 


0.0000 


Jy/pixel 


Look at fit to data - Amplitudes 


eh.plotall obs im compare(obs, out, "uvdist", "amp") 





uvdist lelO 


Look at fit to data — Closure Phase 


eh.plot cphase obs im compare(obs, out, "ALMA","SMA","LMT") 


ALMA - SMA - LMT 


Phase (deg) 


Closure 





6 
GMT (h) 


Step 7: Participate in the Imaging 
Challenge! 


New challenge out now! 
Deadline: December 9 , 2016 


vlbiimaging.csail.mit.edu/imagingchallenge 


* Blind data that you download 
(in uvfits, oifits, and text files) 


Testing Data and Submission Instructions 


e Sample data with truth images 
to help verify your algorithms 
are working 





* Code to help you get started Data Parameters and Noise Properties 


Challenge Number Source Location Telescopes Total Fiux (Janskys) Noise Property 


